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Abstract 

We consider the problem of projecting a probability measure tt on a 
set At iv of Radon measures. The projection is defined as a solution of the 
following variational problem: 

ird |!/i*(p-7r)||i, 

where h £ is a kernel, SI C and * denotes the convolution opera¬ 

tor. To motivate and illustrate our study, we show that this problem arises 
naturally in various practical image rendering problems such as stippling 
(representing an image with N dots) or continuous line drawing (repre¬ 
senting an image with a continuous line). We provide a necessary and suf¬ 
ficient condition on the sequence (7V(jv)iV6N that ensures weak convergence 
of the projections (pJ^jjvgN to tt. We then provide a numerical algorithm 
to solve a discretized version of the problem and show several illustrations 
related to computer-assisted synthesis of artistic paintings/drawings. 


1 Introduction 

Digital Halftoning consists of representing a grayscale image with only black 
and white tones [80) . For example, a grayscale image can be approximated 
by a variable distribution of black dots with over a white background. This 
technique, called stippling, is the cornerstone of most printing digital inkjet 
devices. A stippling result is displayed in Figure [^. The lion in Figure 
can be recognized from the dotted image shown in Figure]^. This is somehow 
surprising since the differences between the pixel values of the two images are far 
from fzero. One way to explain this phenomenon is to invoke the multiresolution 
feature of the human visual system [71121]. Figures an d[TJi are blurred versions 

of Figures and[^ respectively. These blurred images correspond to low-pass 
versions of the original ones and are nearly impossible to distinguish. 

Assuming that the dots correspond to Dirac masses, this experiment suggests 
placing the dots at locations pi,... ,pjv corresponding to the minimizer of the 
following variational problem: 

min 
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Figure 1: Explanation of the stippling phenomenon. Images (a) and (b) are 
similar while the norm of their difference is large. Figures (c) and (d) are 
obtained by convolving (a) and (b) with a Gaussian of variance equal to 3 
pixels. After convolution, the images cannot be distinguished. 


where C denotes the image domain, 5p. denotes the Dirac measure at 

point Pi S TT denotes the target probability measure (the lion) and h is 

a convolution kernel that should depend on the point spread function of the 
human visual system. By letting 

denote the set of iV-point measures, problem Q rereads as a projection problem: 

min ||/i* (tt — u)||? . (3) 

This variational problem is a prototypical example that motivates our study. As 
explained later, it is intimately related to recent works on image halftoning by 
means of attraction-repulsion potentials proposed in [MlEHlIls]. In references 
[nmnn] this principle is shown to have far reaching applications ranging from 
numerical integration, quantum physics, economics (optimal location of service 
centers) or biology (optimal population distributions). 
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In this paper, we extend this variational problem by replacing with 

an arbitrary set of measures denoted A4n- In other words, we want to approx¬ 
imate a given measure tt by another measure in the set Mn- We develop an 
algorithm to perform this projection in a general setting. 

To motivate this extension, we consider a practical problem; how to perform 
continuous line drawing with a computer? Continuous line drawing is a starting 
course in all art cursus. It consists of drawing a picture without ever lifting the 
paintbrush from the page. Figurej^shows two drawings obtained with this tech¬ 
nique. Apart from teaching, it is used in marketing, quilting designs, steel wire 
sculptures, connect the dot puzzles,... A few algorithms were already proposed 
in HOI [Ml HU 0131] • We propose an original solution which consists of setting 
AIat as a space of pushforward measures associated with sets of parameterized 
curves. 

Apart from the two rendering applications discussed in this paper, this paper 
has potential for diverse applications in fields such as imaging, finance, biology,... 



Figure 2: Two examples of continuous line drawing, (a) A sketch of Marylin 
Monroe by Pierre Emmanuel Godet http: //pagazine. com/ using a continuous 
line. A close inspection reveals that the line represents objects and characters, 
(b) Meisje met de Parel, Vermeer 1665, represented using a spiral with variable 
width. Realized by Chan Hwee Chong http;//www.behance.net/Hweechong 

The remaining of this paper is structured as follows. We hrst describe the 
notation and some preliminary remarks in Section We propose a mathemati¬ 
cal analysis of the problem for generic sequences of measures spaces {M.n)n£'n 
in Section]^ In particular, we give conditions on h ensuring that the mapping 
/r I—||/i*/x ||2 defines a norm on the space of signed measures and provide neces¬ 
sary and sufficient conditions on the sequence {M.N)NeN ensuring consistency 
of the projection problem. We propose a generic numerical algorithm in Sec¬ 
tion]^ and derive some of its theoretical guarantees. In Section we study the 
particular problem of continuous line drawing from a mathematical perspective. 
Finally, we present some results in image rendering problems in Section]^ 
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2 Notation and preliminaries 

In this paper, we work on the measurable space (fl, E), where denotes 

the torus An extension to other spaces such as or [0, l]'^ is 

feasible but requires slight adaptations. Since drawing on a donut is impractical, 
we will set = [0,1]'^ in the numerical experiments. 

The space of continuous functions on is denoted C(n). The Sobolev space 
where to S N, is the Banach space of d dimensional curves in 
Q with derivatives up to the m-th order in ^^([0, T]). Let Ma denote the space 
of probability measures on O, i.e. the space of nonnegative Radon measures p 
on fl such that p{fl) = 1. Throughout the paper tt G A4a will denote a target 
measure. Let Ai denote the space of signed measures on with bounded total 
variation, that is p, = p+ — p- where /i+ and p- are two finite nonnegative 
Radon measures and ||/i||TV = M+(^) + M-(^) < oo- 

Let h : n ^ M. denote a continuous function. Let p G A4 denote an arbitrary 
finite signed measure. The convolution product between h and p is defined for 
all a; S n by: 


p-kh{x):= / h{x-y)dp{y) (4) 

Jn 

= p{h{x - •)) 

In the Fourier space, the convolution Q translates to, for all ^ G Z'^ (see e.g., 

my- 

p-kh{^) = pioHo, 


where p is the Fourier-Stieltjes series of p. The Fourier-Stieltjes series coeffi¬ 
cients are defined for all ^ Gl^ by: 


JQ. 


We recall the Parseval formula: 


\h{x)'^ dx = 


Y. Rf) 


Let J : M” —)■ M denote a function and dJ its limiting-subdifferential (or 
simply subdifferential) [HHI]. Let C C M" denote a closed subset. The indicator 
function of C is denoted ic and defined by 

j C^,') — / ® if X G C, 
i-CKX) s otherwise. 


The set of projections of a point xq G K" on C is denoted Pc{xo) and defined 

by 

Pc{xo) = Argmin||x - xq\\1. 

xGC 

The notation Arg min stands for the whole set of minimizers while arg min 
denotes one of the minimizers. Note that Pc is generally a point-to-set mapping 
except if C is convex closed, since the projection on a closed convex set is unique. 
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The normal cone at a: € M" is denoted Nc{x). It is defined as the limiting- 
subdifferential of ic at a;. A critical point of the function J + ic is a point 
X* that satisfies 0 e dJ{x*) + Nc{x*). This condition is necessary (but not 
sufficient) for x* to be a local minimizer of J -I- ic- 


3 Mathematical analysis 


Let 

AfhifJ,) ■■= \\h* (5) 

In this section, we study some basic properties of the following projection prob¬ 
lem: 

min ( 6 ) 

where {M.N)Nen denotes an arbitrary sequence of measures sets in AIa- 


3.1 Norm properties 

We first study the properties of Afh on the space M of signed measures with 
bounded total variation. The following proposition shows that it is well defined 
provided that h e C(n). 

Proposition 1. Let h € C(0) and fxG Ai. Then h-kfiG 

Proof. It suffices to remark that Va; € fl, \h-k fi{x)\ < ||Ai||Tv||/i||oo < +oo. 
Therefore, h-k jj, G L°°{Lt). Since Lt is bounded, h G L°°{Lt) implies that h G 
L^{n). □ 

Remark 1. In fact, the result holds true for weaker hypotheses on h. If h G 
the set of bounded Borel measurable functions, hk fj,G if {LI) since 


Vx e fl, \hk p,{x)\ < IImIItv" 


sup \ h{x)\ 




< - 1 - 00 . 


Note that the L°°-norm is defined with an ess sup while we used a sup in the 
above expression. We stick to h G C{Ll) since this hypothesis is more usual when 
working with Radon measures. 


The following proposition gives a necessary and sufficient condition on h 
ensuring that Afh defines a norm on A4. 

Proposition 2. Let h G C{L1). The mapping Afh defines a norm on A4 if and 
only if all Fourier series coefficients h{f) are nonzero. 

Proof. Let us assume that /i(^) 7 ^ 0, The triangle inequality and 

absolute homogeneity hold trivially. Let us show that ^ 7 ^ 0 ^ Afh{pf) 7 ^ 0. The 
Fourier series of a nonzero signed measure p is nonzero, so that there is ^ G 
such that fi{f) 7 ^ 0. According to our hypothesis h{f) 7 ^ 0, hence p,kh{f) 7 ^ 0 
and Afhipf) ^ 0 . 

On the contrary, if there exists ^0 € such that /i(^o) = 0- The non-zero 
measure defined through its Fourier series by 


A(C) 


1 ife = 6 

0 otherwise 
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satisfies Nh{y) = 0 and belongs to □ 

From now on, owing to Proposition we will systematically assume - some¬ 
times without mentioning - that h G C(fl) and that h(^) ^ 0, G Finally, 
we show that A/^, induces the weak topology on A4. Let us first recall the 
definition of weak convergence. 


Definition 1. A sequence of measures {fJ-N)Ne'M said to weakly converge to 
HGM, if 


lim / f{x)dfiN{x) = / f{x)dfi{x) 

Jn 


for all continuous functions / : —>■ K. The shorthand notation for weak 

convergence is 


hN 


N—^oo 


fl. 


Proposition 3. Assume that h G C{fl) and that h{f) 0, Wf G Then for 
all sequences {nN)N&i in A4 satisfying ||Atjv||Ty < M < -Poo, \/N G N, 


lim AfhihN) = 0 0 . 

N—>-oo N-^-OO 

Proof. Let be a sequence of signed measures in Ai. 

If fj.N 0, then p-NiO = ) —)■ 0 for all f G Z"^. Since |AAf(C)MOI — 

2M\h{^)\ for all f G’L'^ and ^ \2Mh{f)'\^ < oo, dominated convergence yields 

that AfhihN) —>■ 0 . 

Conversely, assume that AfhihN) -G 0. Since the fiN are bounded, there are 
subsequences that converge weakly to a measure v that depends on the 
subsequence. We have to prove that v = Q iov all such subsequences. Since 
AfhihN) — >■ 0, we have pN^f) —>■ 0 for all f G TA. Therefore, v{ff) = 0, G 
This is equivalent to = 0 (see e.g. [161 P-36]), ending the proof. □ 


3.2 Existence of solutions 

The first important question one may ask is whether Problem Q admits a 
solution or not. Theoremprovides sufficient conditions for existence to hold. 

Proposition 4. If At jv is weakly eompact, then Problem § admits at least 
a solution. In particular, if Mn is weakly closed and bounded in TV-norm, 
Problem § admits at least a solution. 

Proof. Assume At at is weakly compact. Consider a minimizing sequence /i„ G 
Mn- By compacity, there is a /r G At at and a subsequence {pLn,f)kGn such 
that pink T- By Proposition [S] Afh induces the weak topology on any 

fe—>^+oo I 

TV-bounded set of signed measures, so that lim A/)i( 7 r — pk) = Afhi'n' — p). 

k—^oo 

Since closed balls in TV-norms are weakly compact, any weakly closed TV- 
bounded set is weakly compact. □ 

A key concept that will appear in the continuous line drawing problem is that 
of pushforward or empirical measure [4] defined hereafter. Let {X, 7 ) denote an 
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arbitrary probability space. Given a function p : X ^ fl, the empirical measure 
associated with p is denoted p* 7 . It is defined for any measurable set B by 

p*j{B) :=7(p-i(B)), 

where 7 denotes the Lebesgue measure on the interval [0,1]. Intuitively, the 
quantity p^j{B) represents the “time” spent by the function p in B. Note 
that p *7 is a probability measure since it is positive and p* 7 (n) = 1. Given a 
measure p. of kind p = p* 7 , the function p is called parameterization of p. 

Let V denote a set of parameterizations p : X ^ Q and M{V) denote the 
associated set of pushforward-measures: 

M{V) := {p = p* 7 ,p e V}. 

In the rest of this paragraph we give sufficient conditions so that a projection 
on M.{V) exists. We first need the following proposition. 

Proposition 5. Let {pn)nen denote a sequence in V that converges to p point- 
wise. Then {pn*j)nen converges weakly to p* 7 . 

Proof. Let / G C{Ll). Since 11 is compact, / is bounded. Hence dominated 
convergence yields f{Pn {x)) - f{p{x))d"i{x) 0. □ 

Proposition 6. Assume that V is compact for the topology of pointwise con¬ 
vergence. Then there exists a minimizer to Problem Q with M.n = M.{V). 

Proof. By Proposition it is enough to show that M.{V) is weakly compact. 
First, A4{V) is bounded in TV-norm since it is a subspace of probability mea¬ 
sures. Consider a sequence in V such that the sequence (Pn* 7 )ra 6 N 

weakly converges to a measure p. Since V is compact for the topology of point- 
wise convergence, there is a subsequence {pnk)keN converging pointwise top GV. 
By Proposition]^ the pushforward-measure p *7 = p so that p G M.{V) and V 
is weakly closed. □ 

3.3 Consistency 

In this paragraph, we consider a sequence (M.N)Ne.n of weakly compact subsets 
of At A- By Proposition 1^ there exists a minimizer p*j^ G At at to Problem ^ 
for every N. We provide a necessary and sufficient condition on (AtArj^veN for 
consistency, i.e. p% tt. In the case of image rendering, it basically means 

N—>-oo 

that if N is taken sufficiently large, the projection p*j^ and the target image tt 
will be indistinguishable from a perceptual point of view. The first result reads 
as follows. 

Theorem 1. The following assertions are equivalent: 
i) For all tt G AIaj P% tt. 

N—^oo 

a) 'is weakly dense in Ma- 

Proof. We first prove ii) => i). Assume that UNenAiN is weakly dense in AIa- 
This implies that, Vtt G At a, ^{pN)NeN G (AtA/jAfeN such that p^ tt. 


7 


From Propositional this is equivalent to lim Afhil^N — t^) =0. Since /i^ is the 

I_I N^oo 

projection 

0 < AfhifJ-N - tt ) < AfhipN - tt ) -)■ 0 . 


Proposition ^ implies that /ilh ^ tt. 

\ I N—^oo 

The proof of i) ^ ii) is straightforward by contraposition. Indeed, if 
UATeN-^TV is not weakly dense in A^a, there exists ttq G Ma that can not 
be approximated weakly by any sequence (/iAr)AreN S CH 


We now turn to the more ambitious goal of assessing the speed of convergence 
of to TT. The most natural metric in our context is the minimized norm 
~ ^)- However, its analysis is easy in the Fourier domain, whereas all 
measures sets in this paper are defined in the space domain. We therefore prefer 
to use another metrization of weak convergence, given by the transportation 
distance. Moreover we will see in Theorem that the transportation distance 
defined below dominates Aft- 

Definition 2. The transportation distance, also known as Kantorovitch or 
Wasserstein distance, between two measures with same TV norm is given by: 

Wi{n,iy) := inf J \\x - y\\.^^dc{x,y) 

where the infimum runs over all couplings of p, and v, that is the measures c 
on n X n with marginals satisfying c(H, H) = p,{A) and c(H,H) = ir^A) for all 
Borelians A. 

Equivalently, we may define the distance through the dual, that is the action 
on Lipschitz functions: 

Wi{fi,iz)= sup yi{f)-ix{f). (7) 

f:Lzp{f)<l 

We define the point-to-set distance as 

Wi{Mn,'!^) ■= inf Wi{p,tt). 

Obviously this distance satisfies: 

Wi{M.N,7r) < Sn := sup inf ITi(/r,7r). (8) 

Theorem 2 . Assume that h G C(i2) denote a Lipschitz continuous function 


with Lipschitz constant L. Then 


Nh{p.-TT) < LWi {qi, tt) 

(9) 

and 


Afh[ki*N - tt) < LWi{Mn,t^) < L5n- 

(10) 


Proof. Let Tx : h{-) i—> h{x — •) denote the symmetrization and shift operator. 
Let us first prove inequality 


where we used the dual definition Q of the Wasserstein distance to obtain the 
last inequality. 

Let/iAT denote a minimizer of inf tt). If no minimizer exists we may 

take an e-solution with arbitrary small e instead. By definition of the projection 
we have: 


- tt) < - tt) < W{hn,t^) < Sn- 


( 11 ) 

□ 


Even though the bound (101 is pessimistic in general, it provides some insight 


on which sequences of measure spaces allow a fast weak convergence. 


3.4 Application to image stippling 

In order to illustrate the proposed theory, we first focus on the case of iV-point 
measures defined in Eq. This setting is the standard one considered 

for probability quantization (see [T^ [TH] for similar results). As mentioned 
earlier, it has many applications including image stippling. Our main results 
read as follows. 

Theorem 3. Let h denote an L-Lipschitz kernel. The set of N-point measures 
satisfies the following inequalities: 


5n 


sup inf IEi(^,7r) < 
ttCAIa 



1 

Afl/d _ 1 


( 12 ) 


and 


sup inf NhitL — 71") < L 
itGA^a 



1 

1 ' 


(13) 


As a direct consequence, we get the following corollary. 

Corollary 1. Let = At(fl^) denote the set of N-point measures. Then 
there exist solutions pLfj to the projection problem §. 

Moreover, for any L-Lipschitz kernel h G C(0); 


a) AfhifJ-N - tt) = O (^LN d'^ . 

Proof. We first evaluate the bound 5n defined in ([^. To this end, for any given 
TT, we construct an explicit sequence of measures piQ,..., the last of which 
is an A^-point measure approximating tt. 

Note that can be thought of as the unit cube [0,1)^^. It may therefore be 
partitioned in smaller cubes of edge length 1/C with C = . We let 

{‘dJi)i<i<cd denote the small cubes and Xi denote their center. We assume that 
the cubes are ordered in such a way that uJi and oJi+i are contiguous. 
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cd 

We define T^{y!i)^xi- The measure /tq satisfies 

i=l 


Wi(7r,^o) ^ - supDiameter(a;i) 

^ i 

/ Vd 1 

^ 2 

but is not an iV-point measure since N'niuJi) is not an integer. 

To obtain an A^-point measure, we recursively build /i/ as follows: 

m{{xi}) = ^ [Nni_i{{xi})\ , 
fii{{xi+i}) = fii-i{{xi+i,xi}) - ^ [Nfii_i{{xi})\ 

iii< {1/cy-i, 

Mixi}) = a {1,1 + 1}. 

We stop the process for I = (l/C)"^ and let jl = fi(^i/c)d. Notice that Nfii{xi) is 
an integer for all i ^ 1 and that /r/ is a probability measure for all 1. Therefore 
fi is an Wpoint measure. Moreover: 

^ ^||a;, - a:;+i||2 
1 

^ iV(7Vi/d_ !)■ 

Since the transportation distance is a distance, we have the triangle inequality. 
Therefore: 


N 




1=1 


Vd 1 

TfVi/d - 1 ' "’Ar(AT/rf-l) 


N- 



The inequality (13) is a direct consequence of this result and Proposition]^ 

We now turn to the proof of Corollary [T] To prove the existence, first notice 
that the projection problem can be recast as 0- Let p = (pi, • • • ,pn) e 

Problem ([^ 




is continuous. 


. The mapping p i—>■ 
therefore consists of minimizing a finite dimensional continuous function over a 
compact set. The existence of a solution follows. Point ii) is a direct consequence 
of Theorem]^ and bound ( [T^ . Point i) is due to the fact that metrizes weak 
convergence, see Proposition]^ □ 
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4 Numerical resolution 


In this section, we propose a generic numerical algorithm to solve the projection 
problem (|^. We first draw a connection with the recent works on electrostatic 
halftoning [Ml 
son’s problem 


in subsection 4.1 


[Ml in subsection 4.2 


We establish a connection with Thom- 
We then recall the algorithm proposed in 
when Mjy is the set of iV-point measures. Finally, we extend this prin¬ 
ciple to arbitrary measures spaces and provide some results on their theoretical 


performance in section 4.4 


4.1 Relationship to electrostatic-halftoning 

In a recent series of papers [23 [28j [TTl [13] , it was suggested to use electrostatic 
principles to perform image halftoning. This technique was shown to produce 
results having a number of nice properties such as few visual artifacts and state- 
of-the-art performance when convolved with a Gaussian filter. Motivated by 
preliminary results in [M], the authors of |M| proposed to choose the N points 
locations p = {pi)i<i<N G as a solution of the following variational problem: 


1 


2iV2 


N N \ ^ r 

EE Hi^Pi Pj) ^ ^ ^ I H{x Pi) diri^x)^ (11) 

i=l 1 = 1 ' i=l 


Repulsion potential 


Attraction potential 


and then extended to 


where H was initially defined as H{x) = ||a :||2 in [M] 
a few other functions in m- The attraction potential tends to attract points 
towards the bright regions of the image (regions where the measure tt has a large 
mass) whereas the repulsion potential can be regarded as a counter-balancing 
term that tends to maximize the distance between all pairs of points. 

Proposition]^ below shows that this attraction-repulsion problem is actually 
equivalent to the projection problem ®> on the set of A^-point measures defined 
in & We let V* denote the set of solutions of ( |I4| and A4(V*) = {p = 
'^p*) P* ^ ’P*}- We also let A4* denote the set of solutions to problem 


Proposition 7. Let h G C(fi) denote a kernel such that |/i|(0 > 0, G Z‘^. 
Define H through its Fourier series by = |/).p(^). Then problems (§ and 
(141 yield the same solutions set: 


M* =M{r*). 


Proof. First, note that since FI and h are continuous both problems are well 
defined and admit at least one solution. Let us first expand the L^-norm in §: 

^||/i*(/r-7r)||^ = i(/i*(/i-7r),/i*(/r-7r)) 

= \{H {p - tt), p - tt) 

= \ {{H*p,p) - 2(iJ*/i,7r) -k (iL *7r,7r)) . 
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Therefore 


Argmin^||/i-*-(/r-7r)||2 = Argmin ^ ((i/- 2(iJ */r, tt)) . 
^ aiGAIjv ^ 


To conclude, it sufhces to remark that for a measure /r of kind M = ;^ Sili ^Pi j 


1 

2 


{{H-kn,fi) - 2{H-kfi,TT)) 


1 


N N 




1 , ^ f 

— ^ / H{x-pi)dTT{x). 

^ i=i 


□ 


Remark 2. It is rather easy to show that a sufficient condition for h to be 
continuous is that H G C^{ri) or H be Holder continuous with exponent a > 2. 
These conditions are however strong and exclude kernels such as H{x) = ||a::|| 2 - 
From Remark^ it is actually sufficient that h G C°°{n) for Mh to be well 
defined. This leads to less stringent conditions on H. We do not discuss this 
possibility further to keep the arguments simple. 


Remark 3. Corollary^ sheds light on the approximation quality of the mini- 
mizers of attraction-repulsion functionals. Let us mention that consistency of 
problem (14) was already studied in the recent papers EH 0123^. To the best 
of our knowledge, Corollary^ is stronger than existing results since it yields a 
convergence rate and holds true under more general assumptions. 

Though formulations ^ and (141 are equivalent, we believe that the pro¬ 
posed one § has some advantages: it is probably more intuitive, shows that the 
convolution kernel h should be chosen depending on physical considerations and 
simplifies some parts of the mathematical analysis such as consistency. However, 
the set of admissible measures A4(H^) has a complex geometry and this formu¬ 
lation as such is hardly amenable to numerical implementation. For instance, 
is not a vector space, since adding two iV-point measures usually leads 
to (2iV)-point measures. On the other hand, the attraction-repulsion formula¬ 
tion ( |14[ ) is an optimization problem of a continuous function over the set . 
It therefore looks easier to handle numerically using non-linear programming 
techniques. This is what we will implement in the next paragraphs following 
previous works |26l EH] ■ 


4.2 Link with Thomson’s problem 

Before going further into the design of a numerical algorithm, let us first show 
that a specific instance of problem § is equivalent to Thomson’s problem [29| . 
This is a longstanding open problem in numerical optimization. It belongs to 
Smale’s list of mathematical questions to solve for the XXIst century [27]. A 
detailed presentation of Thomson’s problem and its extensions is also proposed 
in [Ti] . 

Let E> = {p G ||p ||2 = 1} denote the unit 3-dimensional sphere. Thom¬ 
son’s problem may be enounced as follows: 

Find pG Argmin (15) 

(P1.....P„)GS« ^ \\Pi-Pjh 
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The term J2i^j ||p,_^p.||^ represents the electrostatic potential energy of N elec¬ 
trons. Thomson’s problem therefore consists of finding the minimum energy 
configuration of N electrons on the sphere S. 

To establish the connection between (H and ( [T^ , it suffices to set H{x) = 
n = S and TT = 1 in Eq. (14). By doing so, the attraction potential has 
the same value whatever the points configuration and the repulsion potential 
exactly corresponds to the electrostatic potential. 

This simple remark shows that finding global minimizers looks too ambitious 
in general and we will therefore concentrate on the search of local minimizers 
only. 


4.3 The case of iV-point measures 


In this section, we develop an algorithm specific to the projection on the set of 
A'^-point measures defined in ([^. This algorithm generates stippling results such 
as in Fig.[^ In stippling, the measure is supported by a union of discs, i.e., a sum 
of diracs convoluted with a disc indicator. We simply have to consider the image 
deconvoluted with this disc indicator as tt to include stippling in the framework 
of Wpoint measures. We will generalize this algorithm to arbitrary sets of 
measures in the next section. We assume without further mention that H{^) 
is real and positive for all This implies that H is real and even. Moreover, 
Proposition implies that problems (|^ and (14) yield the same solutions sets. 
We let p = (pi,... ,pn) and set 

N N I ^ r 

H{Pi Pj) ^ E/„ H{x — Pi) dTr{x). (16) 


i=i j=i 


i=l ' 


F(p) 


G(p) 


The projection problem therefore rereads as: 

min J(p). 
pen™ 


(17) 


For practical purposes, the integrals in G{p) first have to be replaced by numer¬ 
ical quadratures. We let G{p) ~ G{p) denote the numerical approximation of 
G{p). This approximation can be written as 

^ N n 

‘^^P^ = iv E E - Pi)'^F 

i=l j=l 

where n is the number of discretization points Xj and Wj are weights that depend 
on the integration rule. In particular, since we want to approximate integration 
with respect to a probability measure, we require that 

71 

E 

i=i 

In our numerical experiments we use the rectangle rule. We may then take tt^ 
as the integral of tt over the corresponding rectangle. After discretization, the 
projection problem therefore rereads as: 


min J(p) := F[p) — G{p). 
pen™ 


(18) 
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The following result [TJ Theorem 5.3] will be useful to design a convergent 
algorithm. We refer to [T] for a comprehensive introduction to the definition of 
Kurdyka-Lojasiewicz functions and to its applications to algorithmic analysis. 
In particular, we recall that semi-algebraic functions are Kurdyka-Lojasiewicz 

m- 

Theorem 4. Let K : K" —>■ M fee function whose gradient is L-Lipschitz 
continuous and let C be a nonempty closed subset o/M". Being given e G (O, 
and a sequence of stepsizes 7 ^^^ such that e < 7 ^^^ < — e, we consider a 

sequence that complies with 


a;(fc+i) e , with G C (19) 


If the function K + ic is a Kurdyka-Lojasiewicz function and if is 

bounded, then the sequence converges to a critical point x* in C. 


A consequence of this important result is the following. 

Corollary 2. Assume that H is a semi-algebraic function with L-Lipschitz 
continuous gradient. Set 0 < 7 < Then the following sequence converges to 


a critical point of problem (18) 


p{k+l) g ^ ^(0) g 


( 20 ) 


If H is convex, 0 < 7 < ^ ensures convergence to a critical point. 

Remark 4. The semi-algebraicity is useful to obtain convergence to a critical 
point. In some cases it might however not be needed. For instance, in the case 
where C is convex and closed, it is straightforward to establish the decrease of 
the cost function assuming only that V J is Lipschitz. Nesterov in 1231 Theorem 
3] also provides a convergence rate in O ^ terms of objective function 

values. 


Proof. First notice that J is semi-algebraic as a finite sum of semi-algebraic 
functions. 

Function J is by Leibniz integral rule. Let dk denote the derivative with 
respect to pk. Then, since H is even 

1 ^ 

dkF{p) = ^ XI “ Pi) (21) 

i=l 

and 

1 " 

dkGijf) = -^ X Wj^H{xj - pk)TTj. (22) 

1=1 
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For any two sets of N points = (p[.^^) 

N 

llVF(pti)) _ Vi^(p(2))||2 = ^ ||a^F(p^i)) - d,F{p^^^) 


fc=l 


2 
2 


. N N ^ 

= EIIE ™(p“-!>")-- p™)" 


fc = l 2 = 1 
iV TV 


< 


pE(Ei|l!>“-p"-(!>?’-!>™)ll»)’ 


< 


iV4 

L' 


fc = l 2=1 
2 N N 




k—1 2=1 

N N 


< 


iV4 


E"(E(ii!>"’-!>f’ii^+ibr’-i'f’ib)’ 


2 L' 


k—1 2=1 

TV N 


-pni+\\pr-p?^\\l 


k—12=1 


=^iip<‘>-p®i& 


and 


||VG(p(i)) - VG(p(2))||2 = V \\dkG{p(^^'>) - dkGip^m 

II Il2 

fc=l 

. N n 2 

= ]V2 ^ II ^ ^ IL 

k^l i=l 

^ N n 2 

<i^E(E>"m-^iid'>-pf‘ii) 

k^l i=l 

7-2 

= - P^^^\\2 

i=i 

Finally, 

||VJ(p(i))-VJ(p(2))||2 

|lVi^(p(i)) _ vF(p(2))|j2 + ||VG(p(i)) - VG(p(2))||2 

<(f + ^)||P<‘>-P®l |2 = fl|p«-P™l| 2 . 

Now, if we assume that H is convex and (this hypothesis is not necessary, 
but simplifies the proof). Then F and G are also convex and C^. We let 
denote the Hessian matrix of F. Given the previous inequalities, we have 
0 ^ V^F ^ ^Id and 0 ^ V^G ^ ^Id- Hence, the largest eigenvalue in 
magnitude of V^(F — G) is bounded above by 

Moreover, the sequence (x^^^)fegN is bounded since is bounded. □ 
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4.4 A generic projection algorithm 

We now turn to the problem of finding a solution of ([^, where At n denotes our 
arbitrary measures set. In the previous paragraph, it was shown that critical 
points of J + iQN could be obtained with a simple projected gradient algorithm 
under mild assumtpions. Although this algorithm only yields critical points, 
they usually correspond to point configurations that are visually pleasing after 
only a few hundreds of iterations. For instance, the lion in Figure was 
obtained after 200 iterations. Motivated by this appealing numerical behavior, 
we propose to extend this algorithm to the following abstract construction: 


1 . Approximate AI_/v by a subset An of n-point measures. 


2 . 


Use the generic Algorithm (19) to obtain an approximate projection fi 
on An- 


* 

n 


3. When possible, reconstruct an approximation /rjv G AIjv of a projection 
using 

To formalize the approximation step, we need the definition of Hausdorff 
distance: 

Definitions. The Hausdorjf distance between two subsets X andY of a metric 
space (M, d) is: 

'Hd{X,Y) :=max< sup inf d{x,y), sup inf d{y,x) > . 

yeY^ex j 

In words, two sets are close if any point in one set is close to at least a 
point in the other set. In this paper, the relevant metric space is the space of 
signed measures A4 with the norm JVh- The corresponding Hausdorff distance 
is denoted . 

The following proposition clarifies why controlling the Hausdorff distance is 
relevant to design approximation sets An- 

Proposition 8 . Let An and AIat be two TV-bounded weakly closed sets of 
measures such that 'Hjy^{An,AiN) < £• Let /i* be a projection on An- Then 
there is a point p-N G Mn such that JVh{yn ~ Tn) < e and Afh{TT — yx) < 
inf A/ji(7r —/r) + 2e. 

N 

Corollary 3. If lim A4n) = 0, then {fj,^)neN converges weakly along 

n—>-oo 

a subsequence to a solution plfi of Problem Q . 

Proof. We first prove Proposition]^ Since An and Mn are bounded weakly 
closed, by Proposition]^ there exists at least one projection /r* on An and one 
projection on AIat. 

Moreover since An and Xi n are bounded weakly closed, they are also closed 
for Mh, so that the infimum in the Hausdorff distances are attained. Hence there 
exists pin G An such that Afhihn — hx) — T~LXh{-^n,Xix) < £ and p-x G Mx 
such that Afhipx—Pn) ^ The proposition follows from the triangle inequality: 

J^hihx - Tt) < Nh[px - Tn)+^h{Pn “ 7^) 

< £+Uh{Pn - tt) 

< £ + Afh{pn ~ Mjv) Y Afhip’x ~ tt) 

^ -^h{p% — tt) + 2£. 
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For the corollary, let us consider the sequence (/r* )„gN as n tends to infinity. 
Since all are in A^a, which is weakly compact, we have a subsequence that 
converges to Since Aft is a metrization of weak convergence on AAn, this 
fj,'^ is indeed a solution to Problem 

- tt) = lim AfhifJ.1, - tt) 

n—^oo 

= inf Mh{TT-^). 


To conclude this section, we show that it is always possible to construct an 
approximation set An Q Al(r2”) with a control on the Hausdorff distance to 
Mn- Let denote an e-enlargement of Aljv w.r.t. the A/^-norm, i.e.: 

Af^ = — I-In) ^ e}- (23) 

We may define an approximation set as follows: 

Al=M{n^)r\M%. (24) 

For sufficient large n, this set is non-empty and can be rewritten as 

= |m = ^ withp= (p,)i<i<„ e, (25) 

where the parameterization set Vn depends on Afjv and e. With this dis¬ 
cretization of Af A at hand, one can then apply (at least formally) the following 
projected gradient descent algorithm: 

pik+i) g p.p, (p{k) _ ^ with e r^. (26) 

The following proposition summarizes the main approximation result: 

Proposition 9. Assume that h is L-Lipschitz. Set e = -I- n^/d_i o.'nd 

An = An, then 

= o . 

Proof. By construction. An satisfies 


sup inf NhitJ^n - Miv) < e. 

Un&An UN^Mn 


Let be an arbitrary measure in At at. By inequality (12), there exists 
Hn € At(fl") such that Nh{pin — fJ-N) < £■ Therefore fin also belongs to A^- 
This shows that 


sup inf Afhifin - (J-n) < £• 

UN&Mn UneAn 


□ 


The approximation process proposed (241 is non-constructive in the does not 
induce any explicit formula for V^- Moreover, can be an arbitrary set and 
the projection on might not be implementable. We will provide constructive 
approximations for specific measures spaces in Section 
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5 Application to continuous line drawing 

In this section, we concentrate on the continuous line drawing problem described 
in the introduction. We first construct a set of admissible measures Mr that 
is a natural representative of artistic continuous line drawings. The index T 
represents the time spent to draw the picture. We then show that using this 
set in problem ensures existence of a solution and weak convergence of the 
minimizers to any tt € J\4a- We hnish by designing a numerical algorithm 
to solve the problem and analyze its theoretical guarantees. 

5.1 Problem formalization 

Let us assume that an artist draws a picture with a pencil. The trajectory of 
the pencil tip can be defined as a parameterized curve p : [0, T] —)■ ft. The body, 
elbow, arm and hand are subject to non-trivial constraints |21j . The curve p 
should therefore belong to some admissible parameterized curves set denoted 
Vt- In this paper, we simply assume that Vt contains curves with bounded 
first and second order derivatives in L'?([0,T]). More precisely, we consider the 
following sets of admissible curves: 

1. Curves with bounded speed: 

= {pG (Wi’~([0,T]))^ p([0,T]) C L!,||p||oo < ai}, 
where q.\ is a positive real. 

2. Curves with bounded hrst and second-order derivatives: 

= {pG (W2’~([0,r]))^ p([0,T]) C fi,||p|U < ai, 

IIpIIoo < a2|, 

where ai and a 2 are positive reals. This set models rather accurately 
kinematic constraints that are met in vehicles. It is obviously a rough 
approximation of arm constraints. 

3. The proposed theory and algorithm apply to a more general setting. For 
instance they cover the case of curves with derivatives up to an arbitrary 
order bounded in L'^ with g € [l,oo]. We let 

= {p e (W'"’«([0,T]))^ p([0,T]) c fl, 

Vi e {!,...,to}, < a*|. 

where (0^)1=!...m are positive reals. This case will be treated only in 
the numerical experiments to illustrate the variety of results that can be 
obtained in applications. 

Note that all above mentionned sets are convex. The convexity property will 
help deriving efficient numerical procedures. 

In the rest of this section, we consider the following projection problem: 

inf Ab(/r-7r), (27) 
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with a special emphasis on the set A4 since it best describes standard 

kinematic constraints. This problem basically consists of finding the “best” way 
to represent a picture in a given amount of time T. 


5.2 Existence and consistency 

We first provide existence results using the results derived in Section|^for q = oo. 


Theorem 5. For any m € N*, Problem (27) admits at least one solution in 


Proof. From Proposition it suffices to show that is compact for the 

topology of pointwise convergence. 

Let (pn)nGN be a sequence in that converges pointwise to p. Since 

is in its (m — l)-th derivative is Lipschitz continuous. By definition of 

7^™’°°, the pi™ are both uniformly bounded by am-i and Om-Lipschitz, hence 
equicontinuous. Next, by Ascoli’s theorem, up to taking a subsequence, pn™ 
uniformly converges to a continuous Integrating yields that pn'^ —>■ p*^®^ 

uniformly for all i < m — 1, so that ^ cti for i < to — 1. Finally, a 

limit of L-Lipschitz functions is also L-Lipschitz, so that ||p^"*^ || < am- Hence 

p S V^'°°, ending the proof. □ 


Let us now turn to weak convergence. 


Theorem 6. Let T be an arbitrary positive real. Let pip S M (7^™’°°) 
any solution of Problem (27). Then, for any Lipschitz kernel h G C(Q): 


,) ^ ,r. 


iij = O ^, 

Proof. Let us consider a function u : [0,1] —)■ K such that: 

• The TO-th derivative is bounded by Um, that is < Om- 

• For all integers i € {!,...,to — 1}, endpoint values are zero, that is 
mW(o) = u«(l) = 0. 

• Start point is zero, that is m(0) = 0. 

• Endpoint is positive, that is m(1) = C > 0. 


Let X and y in H, such that ||a; — p||2 = Cr™, and let r^y be the unit vector 
from X to y. Then, for r small enough, the function s[a:,p] : t x + Txyu{^) 
belongs to 7^™’°°, with all its first (to — 1) derivatives zero at its endpoints. The 
condition r small enough is for controlling the norm of the i-th. derivatives for 
i < TO — 1, which scale as r™”*. 

Now, let us split H = [0, in N‘^ small cubes uji. We may order them such 
that each uii is adjacent to the next cube Wi+i. We write xi for the center of uji. 
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We now build functions s € V^'°° by concatenating paths from Xi to Xi+i and 
waiting times in xf. 


1 

t? - = f- 





s{t) 


Xi 


if tj <t < tf, 
if 


under the condition T > Tjq := {N^ — 1) (j^) , that is to say that we have 

enough time to loop through all the cube centers. 

Let now tt G Ma- We may choose — tj < TTr{uJi) for all i. Then, we may 
couple TT and s*7t with c{xi,uji) = ‘.p * . Since the small cubes have radius 
'/d/N and the big one has radius \/d, we obtain: 


Wi(7r, s*7t) < 


I 

2N ^ T T 

I 


VdT-TN 

m T 




In particular, taking N = , we find that ITi (M {V^’°°) ,tt) = O {t rr,{d+i)-i 

hence IJ-p is weakly dense in Ma- 

□ 


5.3 Numerical resolution 


We now turn to the numerical resolution of problem (27). We first discretize 
the problem. We set At := ^ and define discrete curves s as vectors of 
We let s{i) G denote the curve location at discrete time i, corresponding to 
the continuous time iAt. 

We define £>1 : ^ 

follows: 


the discrete first order derivative operator, as 


(n \C\ ® if f = 1, 

_ 1) iiiG{2,..., N}. 

In what follows, Di denotes a discretization of the derivative operator of order 
i. In the numerical experiments, we set D 2 = —D*Di. 

We define P™’®, a discretized version of as follows: 

= {s G such that Vi G {!,..■ N}, s{i) G fl, (28) 

and Vj G {I,... ,m}, ||Pjs||, < a^}. (29) 

( N-d \ I 

for q G [1,+cxo) and ||a;||oo = 

max ||xj||2. 
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The measures set can be approximated by the set of fV-point 


measures 


^ ). From Corollary]^ it suffices to control the Hausdorff dis¬ 
tance Pvyj(Af to ensure that the solution of the discrete 


problem § with AfTV — *^) ^ §ood approximation of problem ( |2T[ ) . 

Unfortunately, the control of this distance is rather technical and falls beyond 
the scope of this paper for general m and q. In the following proposition, we 
therefore limit ourselves to the case m = l,q = oo. 


Proposition 10. < oi^- 

Proof. 1. Let us show that 


sup 


inf 


) ^ 


Let /i e and denote by p S 


P. 


N-l 


that fj, = p* 7 . Define p = 77 




i=0 


a parameterization such 
y Then a parameterization of 


fi is defined by s(i) = Moreover, for i G {2,... A}, |(Dis)(z)| = 


1 

At 


iT 

A 


- 1 )^ 
A 


1 

At 


Oi- Therefore s G Pb°° ■ 


’ (»-l)T 


p{t) dt 


< 


1 fTT 

At IiPzllL 


\p{t)\ dt < 


Let us consider the transportation map coupling the curve arcs between 
times {i — 1)^ and and the Diracs at p (i^). Then IUi(p* 7 , 5 * 7 ) ^ 

^ 1 

F- 

^ N 


sup 


s(t)-s( 


T 


2. Let us fix /r G Af and let s G such that 3*7 = /i. We set 


p( 0 ) = s(l), and: 
p{t) = 


s(l) for t G 

+ “ *) + l)-s(*)J'or t G 


o>f]. 

iT (i+l)T 


N ’ N 

,tG A-1}. 


Since s G and D is convex, p([0,T]) C D. Moreover, p is continu¬ 
ous and piecewise differentiable. Finally, for i G {1,..., A — 1} and t G 


iT (i+l)T 
N ’ N 


’ p(^) = St ('5(* + 1) - s(*)) = Di{s){i). Therefore, ||p|| 


€ 


Oi, ensuring that p G With the same coupling as above, we have 

VFi(p» 7 , s* 7 ) ^ oi^, which ends the proof. 

□ 


To end up, let us describe precisely a solver for the following variational 
problem: 

inf Afh{fJ.-Tr). (30) 

mgSi(p).-”) 

We let A4* denote the set of minimizers and P* denote the associated set of 
parameterizations. 
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Algorithm 1: A projection algorithm on At . 

Input: 

- tt: target measure. 

- N: a number of discretization points. 

- G initial parameterized curve. 

- H: a semi-algebraic function with Lipschitz continuous gradient. 

- nit: number of iterations. 

Output: 

- an approximation of a curve in V* . 

- an approximation of an element of At*, 
for 0 < A: < nit do 

- Compute 

- Set = Ppcoo 


Remark 5. The implementation of Algorithm^ requires computing the gradi¬ 


ents (21) and (22) and computing a projection on Both problems are 

actually non trivial. 

The naive approach to compute the gradient of F consists of using the explicit 
formula This approach is feasible only for a small amount of points N 

(less than 1000^ since its complexity is O {N"^) ■ In our numerical experiments, 
we therefore resort to fast summation algorithms ]2IA commonly used in 
particles simulation. This part of the numerical analysis is described in \2f^ 
and we do not discuss it in this paper. 

The set and more generally the sets P^’‘‘ are convex for q G [l,oo]. 

Projections can be computed using first-order iterative algorithms for convex 
functions. In our numerical experiments, we use accelerated proximal gradient 
descents on the dual problem precise description is given in m- 


6 Results 

To illustrate the results, we focus on the continuous line drawing problem dis¬ 
cussed throughout the paper. It is performed using Algorithm [T] In the follow¬ 
ing experiments, we set H as a smoothed L^-norm. This is similar to what was 
proposed in the original halftoning papers in [511HH] ■ 
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6.1 Projection onto 

In this part, we limit ourselves to the projection onto as studied in the 

previous section. In Figure we show the evolution of the curve across 
iterations, for different choices of After 30,000 iterations, the evolution 
seems to be stabilized. The cost function during the 400 first iterations is 
depicted in Figure]^ for the three different initializations. 

In Figure we show the projection of the famous Girl with a Pearl Earring 
painting, after 10, 000 iterations. To really see the precision of the algorithm, 
we advise the reader to blink the eyes or to take a printed version of the paper 
away. From a close distance, the curves or points are visible. From a long 
distance, only the painting appears. 

6.2 Projection onto 

We now consider projections onto more general measure spaces, such as M 
in order to show that different measures spaces can be considered. In Fig.|^ we 
show different behaviours for different m £ {1,2} and g £ (1,2, oo}. We also 
show a large scale example with a picture of Marylin Monroe in Figure 


7 Conclusion 

We analyzed the basic properties of a variational problem to project a target 
Radon measure tt on arbitrary measures sets At tv- We then proposed a numer¬ 
ical algorithm to find approximate solutions of this problem and gave several 
guarantees. An important application covered by this algorithm is the projec¬ 
tion on the set of 7V-point measures, which is often called quantization and 
appears in many different areas such as finance, imaging, biology,... To the best 
of our knowledge, the extension to arbitrary measures set is new, and opens 
many interesting application perspectives. As examples in imaging, let us men¬ 
tion open topics such as the detection of singularities [5] (e.g. curves in 3D 
images) and sparse spike deconvolution in dimension d [S] . 

To finish, let us mention an important open question. We provided necessary 
and sufficient conditions on the sequence (AlAr)TVGN for the sequence of global 
minimizers fo weakly converge to tt. In practice, finding the global 

minimizer is impossible and we can only expect finding critical points. One may 
therefore wonder whether all sequences of critical points weakly converge to tt. 
An interesting perspective to answer this question is the use of mean-filed limits 

i- 
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Figure 3: Projection of the lion image onto with N = ^000. The figure 
depicts with several values of the iterate k in Algorithm 1 
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Figure 4: Decay of the cost function J for the three experiments depicted in 
Fig. § We represent log]^Q( J(fc) — m) for k < 400 where m is the mimimal value 
of J during the first 30,000 iterations. 
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Figure 5: Projection of Meisje met de Parel, Vermeer 1665, onto with 

N = 150,000. The figure depicts gfiohoo) obtained with Algorithmj^ 
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Figure 6: Projection of the lion image onto with N = 8, 000, and m € 
{1,2} and q € {1, 2, ooj. 
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Figure 7: Projection of Marylin image, onto the set: 

C = {p € (VF^’°°([0,r]))2, sup (||I?ip(*)|| 2 ) < ai, sup i\\D 2 p{i)\\ 2 ) < 02 }, 
with N = 100,000. The figure depicts s(io>ooo) obtained with Algorithm 
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